######
# to do the analysis of the Montecarlo estimations
# Created: 20161123
# Last modified: 20190830
######
# Calling packages
foo <- function(x){
  for( i in x ){
    #  require returns TRUE invisibly if it was able to load package
    if( ! require( i , character.only = TRUE ) ){
      #  If package was not able to be loaded then re-install
      install.packages( i , dependencies = TRUE )
      #  Load package after installing
      require( i , character.only = TRUE )
    }
  }
}

#  Then try/install packages...
foo(c("plyr","doBy","data.table"))
# Definig Globals 
Root <- "C:/Users/ja.guerra"
#setwd("C:/Users/uctpjag/Dropbox/Temporal")
WD <- paste(Root,"/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes",sep="")
Codes <- paste(WD,"/C_ANALYSIS",sep="")
Results <- paste(WD, "/Results/ReStat/",sep="")
Resultsks <- paste(Results, "/KS/",sep="")
Graphs<- paste(Results,"Graphs/",sep="")
Graphsks <- paste(Graphs, "KS/",sep="")
algo <- c("","KS")
setwd(Results)

i0 <- 2
while(i0 <= 2) {
  setwd(paste(Results,algo[i0],"/",sep=""))
  sink("SummaryMontecarlos.txt")
  
  #JeJ
  #######
  #calling all files that contain Je
  RDatalist <- list.files(path=paste(Results,algo[i0],sep=""),pattern="dataJ+")
  sumJs <- vector(mode="list",length=length(RDatalist))
  summseJs <- sumJs
  for (Np in RDatalist) {
    i <- which(RDatalist==Np)
    load(paste(paste(Results,algo[i0],"/",sep=""),Np,sep=""))
    if(i0 ==2) dataJeJ <- dataJeJks
    numb <- paste(regmatches(Np,gregexpr("[0-9]",Np))[[1]],collapse="")
    print(numb)
    dataJeJ$Np <- Np
    sumJ <- summaryBy(. ~ Np, data=dataJeJ, FUN=c(median,mean,sd), na.rm=TRUE, keep.names=T,order=TRUE)
    sumqJ <- summaryBy(. ~ Np, data=dataJeJ, FUN=c(quantile), probs= c(0.25,0.75), na.rm=TRUE, keep.names=T,order=TRUE)
    sumJ <- cbind(sumJ, sumqJ)
    sumJ <- sumJ[,order(names(sumJ))]
    sumJ <- sumJ[,min(which(colnames(sumJ) %in% grep("expw", names(sumJ),value=TRUE))):length(sumJ)]
    print(sumJ)
    sumJs[[i]] <- sumJ
    true <- dataJeJ[(colnames(dataJeJ) %in% paste("J",c(1:5),sep=""))]
    datamseJ <- cbind(dataJeJ[(names(dataJeJ)  %in% paste("Jexpw",c(1:5),sep=""))]-true,
                      dataJeJ[(names(dataJeJ)  %in% paste("Jsw",c(1:5),sep=""))]-true,
                      dataJeJ[(names(dataJeJ)  %in% paste("Jswk",c(1:5),sep=""))]-true)
    datamseJ <- as.data.frame(datamseJ^2)
    datamseJ$Np <- Np
    summseJ <- summaryBy(. ~ Np, data=datamseJ, FUN=c(mean), na.rm=TRUE, keep.names=T)
    print(summseJ)
    summseJs[[i]] <- summseJ
  }
  sumJs <- rbindlist(sumJs)
  write.csv(sumJs, file = "sumJsinter.csv")
  summseJs <- rbindlist(summseJs)
  write.csv(summseJs, file = "summseJsinter.csv")
  
  #ded
  #######
  #calling all files that contain de
  RDatalist <- list.files(path=paste(Results,algo[i0],sep=""),pattern="datad+")
  sumds <- vector(mode="list",length=length(RDatalist))
  summseds <- sumds
  for (Np in RDatalist) {
    i <- which(RDatalist==Np)
    load(paste(paste(Results,algo[i0],"/",sep=""),Np,sep=""))
    if(i0 ==2) dataded <- datadedks
    numb <- paste(regmatches(Np,gregexpr("[0-9]",Np))[[1]],collapse="")
    print(numb)
    dataded$Np <- Np
    sumd <- summaryBy(. ~ Np, data=dataded, FUN=c(median,mean,sd), na.rm=TRUE, keep.names=T,order=TRUE)
    sumd <- sumd[,order(names(sumd))]
    sumd <- sumd[,min(which(colnames(sumd) %in% grep("expw", names(sumd),value=TRUE))):length(sumd)]
    print(sumd)
    sumds[[i]] <- sumd
    true <- dataded[(colnames(dataded) %in% paste("d",c(2:5),sep=""))]
    datamsed <- cbind(dataded[(names(dataded)  %in% paste("dexpw",c(2:5),sep=""))]-true,
                      dataded[(names(dataded)  %in% paste("dsw",c(2:5),sep=""))]-true,
                      dataded[(names(dataded)  %in% paste("dswk",c(2:5),sep=""))]-true)
    datamsed <- as.data.frame(datamsed^2)
    datamsed$Np <- Np
    summsed <- summaryBy(. ~ Np, data=datamsed, FUN=c(mean), na.rm=TRUE, keep.names=T)
    print(summsed)
    summseds[[i]] <- summsed
  }
  sumds <- rbindlist(sumds)
  write.csv(sumds, file = "sumdsinter.csv")
  summseds <- rbindlist(summseds)
  write.csv(summseds, file = "summsedsinter.csv")
  
  #aea
  #######
  #calling all files that contain ae
  RDatalist <- list.files(path=paste(Results,algo[i0],sep=""),pattern="dataae+")
  sumas <- vector(mode="list",length=length(RDatalist))
  summseas <- sumas
  for (Np in RDatalist) {
    i <- which(RDatalist==Np)
    load(paste(paste(Results,algo[i0],"/",sep=""),Np,sep=""))
    if(i0 ==2) dataaea <- dataaeaks
    numb <- paste(regmatches(Np,gregexpr("[0-9]",Np))[[1]],collapse="")
    print(numb)
    dataaea$Np <- Np
    suma <- summaryBy(. ~ Np, data=dataaea, FUN=c(median,mean,sd), na.rm=TRUE, keep.names=T,order=TRUE)
    suma <- suma[,order(names(suma))]
    suma <- suma[,min(which(colnames(suma) %in% grep("expw", names(suma),value=TRUE))):length(suma)]
    print(suma)
    sumas[[i]] <- suma
    true <- dataaea[(colnames(dataaea) %in% paste("a",c(2:5),sep=""))]
    datamsea <- cbind(dataaea[(names(dataaea)  %in% paste("aexpw",c(2:5),sep=""))]-true,
                      dataaea[(names(dataaea)  %in% paste("asw",c(2:5),sep=""))]-true,
                      dataaea[(names(dataaea)  %in% paste("aswk",c(2:5),sep=""))]-true)
    datamsea <- as.data.frame(datamsea^2)
    datamsea$Np <- Np
    summsea <- summaryBy(. ~ Np, data=datamsea, FUN=c(mean), na.rm=TRUE, keep.names=T)
    print(summsea)
    summseas[[i]] <- summsea
  }
  sumas <- rbindlist(sumas)
  write.csv(sumas, file = "sumasinter.csv")
  summseas <- rbindlist(summseas)
  write.csv(summseas, file = "summseasinter.csv")
  
  i0 <- i0+1
  sink()
}
